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; ABSTRACT 

We present accurate predictions of the correlation function of hotspots in the mi- 
$_i ' crowave background radiation for gaussian theories such as those predicted in most 

inflation models. The correlation function of peaks above a certain threshold depends 
' only on the threshold and the power spectrum of temperature fluctuations. Since there 

■ are both potentially observable quantities in a microwave background map, there are 

fN| | no adjustable parameters in the predictions. These correlations should therefore pro- 

vide a powerful test of the gaussian hypothesis, and provide a useful discriminant 
between inflation and topological defect models such as the cosmic string model. The 
correlations have a number of oscillatory features, which should be detectable at high 



> 



' signal-to-noise with future satellite experiments such as MAP and Planck. 
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1 INTRODUCTION 

The microwave background radiation is principally a relic of the early Universe, containing information on the fluctuations 



O ' present at a redshift of about 1000. Detailed knowledge of the fluctuations is extremely valuable as a tool for distinguishing 
structure formation theories. The main alternative to the microwave background for this purpose is the distribution of galaxies 



at the present day, or at moderate redshifts. The microwave background has two great advantages over galaxy redshift surveys. 
Firstly, one needs to make no assumptions about the relationship between galaxy clustering and mass clustering (the issue 
of bias) , and secondly, the microwave background probes the Universe when the fluctuations were extremely small, so linear 
perturbation theory can describe the predicted fluctuations to high accuracy. Against these advantages must be offset the 



5— | ' problems of foreground contamination and the technical difficulties of producing a high signal-to-noise, high-resolution map of 
5^ i large areas of sky. Future satellite experiments planned by ESA and NASA (Bersanelli et al. 1996, Jungman et al. 1996) will 
produce such maps, and these will allow the spectrum of temperature fluctuations to be measured with high accuracy. The 
primary goal of such experiments will be to determine a number of cosmological parameters, such as the density parameter 
Qo, the Hubble constant, and the mix of species of dark matter, all of which, in inflationary theories, influence the detailed 
shape of the temperature power spectrum. In such theories, the temperature map is a random gaussian field, whose properties 
are determined by the power spectrum (or correlation function) alone. It is vitally important to test whether the map is 
indeed gaussian, as if it is not, the proposed parameter estimation method (e.g. Jungman et al. 1996) is invalid. In addition, a 
non-gaussian map of the primary temperature anisotropies would indicate a quite different structure formation scenario from 
the standard inflation model. 

There are several ways to test the gaussian hypothesis, such as the 3-point function (e.g. Hinshaw et al. 1994, Falk, 
Rangarajan & Srednicki 1993, Luo & Schramm 1993, Gangui et al. 1994), the genus and Euler-Poincare statistic (Coles 1989, 
Gott et al. 1990, Luo 1994b, Smoot et al. 1994), the bispectrum (Luo 1994a, Heavens 1998, Ferreira, Magueijo & Gorski 
1998), studies of tensor modes in the CMB (Coulson, Crittenden & Turok 1994), and peak statistics (Bond & Efstathiou 
1987, Kogut et al. 1995, Kogut et al. 1996). For all of these it is possible in principle to make predictions for gaussian fields, 
and each method has its advantages. The bispectrum is attractive because the error analysis can be done without recourse 
to numerical simulation, and can treat the twin problems of missing sky regions and correlated noise (Heavens 1998), albeit 
at large computational expense. In this paper we explore an alternative: the correlation function of peaks in the microwave 
background. This has the advantage of being readily calculable, and also we find that the predicted correlation functions have 
oscillatory features which are not obviously present in the temperature autocorrelation function. In much the same way as 
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the oscillations in the power spectrum allow high precision determinations of cosmological parameters, the peak correlation 
function oscillations should allow a high precision test of the gaussian hypothesis. In the latter case, the correlation function 
of peaks above any given threshold is completely determined once the power spectrum is known; there are no free parameters. 

In this paper, we calculate the 2-point correlation function of peaks in a 2D gaussian field, for any given power spectrum. 
We make no approximations other than to compute the correlation function for small angular separations so we assume the sky 
is flat and translate the sky power spectrum to an analogous 2D spectrum. We follow the peak statistics methods developed 
initially by Peacock & Heavens (1985) and Bardeen et al. (1986). 

The layout of the paper is as follows. In section 2 we describe the method, and in section 3 present the results and 
compare with an approximate correlation function method introduced by Bardeen et al. (1986). In section 4 we discuss the 
results and the possibilities for Planck and MAP. 



2 METHOD 

In this section, we compute the two-point correlation function of local maxima in 2D gaussian random fields. For small 
separations, this calculation should be accurately applicable to the clustering of peaks in the microwave sky (see Bond & 
Efstathiou (1987) for further discussion of this point). The calculation is laid out as follows: after some definitions, we compute 
the 12 x 12 covariance matrix for the 12 variables which are required to specify the two peaks. We then invert this matrix to 
obtain the multivariate gaussian distribution for the variables, then apply the peak constraint and integrate over the remaining 
variables to obtain the correlation function. 



2.1 Peaks in 2D 

We define the temperature fluctuation by <5(x) = T(x)/T — 1, and its 2D Fourier transform by <5k = J d 2 k<5(x) exp(ik • x). 
The fluctuations are specified entirely by the power spectrum, P(k), defined by 

(S k S k ,) = (2tv) 2 P(k)S D (k + k') (1) 

where angle brackets indicate ensemble averages, and 8 D is a 2D Dirac delta function. We will require moments of the power 
spectrum, defining 

^ = (2^/ d2kP(fc)fe2j ' (2) 
and spectral parameters 

7 = «T?/((ro<7a) 0, = \/2— . (3) 

At peaks, the 2D gradient Si = dS/dxi vanishes, and the eigenvalues of the second derivative matrix Sij = d 2 5 / dxidxj are all 
negative. Since Sij is symmetric, we are concerned with 6 independent numbers v = {5, S x , S y , S xx , S xy , S yy } to specify a local 
maximum. For a gaussian field, the probability density function for the 12 variables v = {v' 1 ', v' 2 '}, evaluated at two points 
separated by a distance r, is 

P(V1 ' V2) = (2^||M||^ ° XP H ViM ^') • (4) 

The covariance matrix is Mij = {(vi — Vi)(vj — Vj)}, overlines indicate ensemble averages, and the summation convention is 
assumed. Note that vt = in our case, and is the i,j component of the inverse of M. 



The number density of peaks is a sum of 2D Dirac delta functions, centred on peaks x. p k, q : 
n pk = 2J S D (x - x pk ,q) (5) 

In the vicinity of a peak, we may expand 

S(x) ~ <5(x pfc ) + - ^2 dx _ dx ( x - x pfcM x - x pfc)j (6) 



SO 



d 2 S \ 1 dS 



and the Dirac delta function in ^ is 

<5 c (x — Xpfc) ~ | det Sij\S D (Si). (8) 
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The mean number density is 
(n pfc (x)} = (j det 5 ZJ \S D (Si)) 

llSdSxX dS X y dSyy \ det <5»j I ^(5, Si = 0, Sij) 

and the range of integration is set by the requirement that the eigenvalues of - 
1 



(9) 



-Sij are all positive. These eigenvalues are 



A 



'{SxX + Syy) ± y/ (S XX — Syy) 2 + ASxy 



(10) 



If we introduce the notation 
S 



x 

Y 

Z 



0"2 



25. 



<J2 



then the conditions for a maximum are X > and Y 2 + Z 2 < X 2 . The determinant is 

J^"2 y-2 ^2 

det Sa — ■ 



4a 2 



(11) 



(12) 



To compute the mean number density, we need the covariance matrix for the 6 variables specifying the field and its first and 
second derivatives. The non-zero components are (u 2 ) = (X 2 ) — 1; (vX) = 7; (Y 2 ) = (Z 2 ) — 1/2; (S 2 .) = (8 2 ) = crf/2. The 
matrix can readily be inverted, and the integrations over Y , Z and X performed. The expression for the differential number 
density as a function of peak height, n v k(y) is given in Bond & Efstathiou (1987), equation (A1.9): 

1 2 



n vk (v) 



where 



G(7, x* 



(2tt) 3 / 2 



(13) 




V3-27 



(14) 



The correlation function is obtained similarly. We find the covariance matrix for the 12 variables, invert it to find their 
probability distribution, and integrate subject to the constraints. The cross-correlation function of peaks of height v\ and V2 
is then given by 

W X 't- Y l fV X 2~ Y 2 



(Xl - Y? - Z\) (Xl 



x 1= o Jx 2 =o Jy 1 =-x 1 Jy 2 = -x 2 J Zi_=-yJX'l-Y'l J z 2 =- v /x%-y. 



dX 1 dX2dY 1 dY 2 dZ 1 dZ2 x 



Y4 



Zi)p(u 1 ,X l ,Yi,Zi,5 : 



(i) 



0,Sy 



(1) 



Q,V2,X2,Y 2 ,Z 2 ,S y x 



(2) 



0,^ 



(2) 



0). (15) 



The correlation function for peaks above a certain threshold v is obtained by adding two further integrations over v\ and ^2, 
and replacing the differential number densities n p k(v) in the denominator of ( [li]) by numerically-evaluated integrals n p k(> v). 
The calculation of the covariance matrix and the joint probability distribution is rather technical, and so appears in an 
appendix. 



3 RESULTS 

In the flat-sky approximation, the power spectrum P(k) which we require is related simply to the conventional power spectrum 
of spherical harmonic coefficients Ce = (|a^ m | 2 ), where 8T(9,4>) = **£2i tm CL£ m Yp(9,(j}) (Bond & Efstathiou 1987): 

P(k) = Ct with k ~ t (16) 

We run CMBFAST (Seljak & Zaldarriaga 1996) to generate the power spectrum Ct, and interpolate to obtain a continuous 
power spectrum P(k). We model the beam with a gaussian of FWHM b, so multiply the power spectrum by a gaussian 
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Figure 1. The correlation function of peaks in the microwave background, for a cold dark matter (CDM) model (top left; baryon density 
Q.g = 0.05, &cdm = 0.95), a mixed dark matter model containing hot dark matter (top right; Qg = 0.05, &CDM = 0.8, Qhdm = 0.15), 
and a fiat cosmological constant model (bottom left; f2g = 0.05, OcDM = 0.15, Q\ = 0.8). A Hubble constant Hq = 60 km s — 1 Mpc - 1 
is assumed and the beam is 5 arcminutes FWHM. At bottom right, we show the effect on the cosmological constant model of altering 
the beam size. 5, 10, 15 and 20 arcminute beams are shown. 



exp [—a 2 £(£ + 1)] , with a = b/V8 In 2. For peaks above a threshold, an 8D numerical integration is performed (although one 
could be done analytically), and, despite the integrand consisting of about 3000 lines of code, each point in the correlation 
function can be computed accurately in tens of seconds on a fast workstation. 

In Fig. ^, we show the peak-peak correlation functions for a variety of models, for a 5 arcminute beam and peaks above 
various thresholds. The striking feature of these graphs is the wealth of structure in the correlation functions. The first peak is 
not due to the beam, rather due to the sharp turndown in the power spectrum itself at around £ = 1500, due to the thickness 
of the last scattering surface and/or Silk damping. There may be a correspondence in the peak positions with analogous 
peaks in the power spectrum, but it could also be ringing due to the abrupt cutoff. Note in particular that the temperature 
autocorrelation function (the Fourier transform of the power spectrum) is very smooth for these separations, without any 
oscillatory features. We also show in Fig. ^ the effect of beam-smearing. Correlation functions for beams of 5, 10, 15 and 20 
arcminute are shown. The oscillations are more numerous and pronounced for smaller beam sizes, so it is advantageous to 
work at high resolution. 

Later in this paper we will want to compute the correlations of estimates of the peak-peak correlation function. This 
is difficult analytically, so we have generated sky realisations, of 3072 x 3072 pixels, each 1 arcminute across, located local 
maxima, and computed the correlation function (see Fig. ^). Edge effects are irrelevant as we wrap the box round for pair 
counting. Fig. ^| shows the average of 10 such realisations, with the error on the mean plotted for each bin. The beam is 5 
arcminutes and the spectrum is from a mixed dark matter model, which we have apodized with a cosine bell to remove power 
at long wavelengths. The apodizing affects k < 30, and is introduced to avoid discreteness errors arising at low-fc from our 
continuum approximation used to evaluate the moments of the power spectrum and the A a/ 3 7 , defined in the Appendix. Note 
the good agreement with the theoretical curve. 

This numerical approach allows us to compute the correlation coefficients of the estimators £ of the peak-peak correlation 
function. In Fig. ^ we show, from 100 realisations of a mixed dark matter model, ( (£, — ) (£,■ — £j ) ) / (<7j <Tj ) where fj is the mean 
value and Oi the variance of the estimates of the correlation function at radius rj. These correlations would be required in order 
properly to quantify the significance of any departure from gaussian statistics. Note, however, that the correlation matrix 
is very close to diagonal, so distinct estimators are almost uncorrelated. We have also tested the approximate correlation 
function considered in 3D by Bardeen et al. (1986) and in 2D by Bond & Efstathiou (1987). In this approximation, more 
analytic progress is possible. The approximation ignores all derivatives of the temperature autocorrelation function, and it 
should be accurate for large separations. In the notation of the appendix, this means setting all the Xc,p v = except for the 
normalised autocorrelation function itself, Aooo- We show an example of the approximation, with the accurate solution, in 
Fig. |^. The approximation does quite well, especially at large separations, but fails to reproduce the oscillations. In Fig. || 
we show the expected errors from MAP and Planck for a mixed dark matter model. The data were obtained by repeatedly 
simulating part of the sky on a grid, and rescaling the errors to simulate 75% sky coverage. Details are given in the caption. 
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Figure 2. A 10° X 10° patch of a mixed dark matter simulated sky, at 5 arcminute resolution characteristic of the highest frequency 
Planck channels. Peaks above la are marked with circles. 




Figure 3. The v = 1 threshold from the mixed dark matter model in Fig. with the correlation function from simulated gaussian 
fields superimposed. Also shown is the approximate correlation function computed by Bond & Efstathiou (1987) by ignoring derivatives 
of the temperature autocorrelation function. 



4 DISCUSSION 

In this paper we have shown that the peak-peak correlation function for the microwave sky can be computed to high accuracy 
for gaussian fields. The only approximation we make is that the sky is flat, but since the features in the correlation function 
appear on scales less than about 2 degrees, the computations should be very accurate. With upcoming high-resolution CMB 
experiments MAP and Planck, these results should allow a high-precision test of the hypothesis that structure grows from 
gaussian initial fluctuations, such as predicted in most inflationary models. We have not computed the expected peak-peak 
correlation function from defect models, but it is clear from the abundance of structure in the correlation function that 
agreement with these predictions, for all peak heights, would constitute a powerful argument in favour of gaussian models. We 
have also investigated numerically the errors expected from MAP and Planck; many of the features in the correlation function 
should be detectable with high significance. We note the similarity in the number of features in the peak-peak correlation 
function as in the power spectrum (although there may not be a simple correspondence between the two). The use of the two 
curves will, however, be rather different. To fit the power spectrum, one has the freedom to adjust fifteen or more parameters, 
under the gaussian hypothesis. With the peak-peak correlation function, there are no free parameters; for a gaussian process, 
all the statistics are specified once the power spectrum is fixed. Thus there is no freedom. Note also that the information 
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Covariance matrix 




10 100 
8/ arcmin 



Figure 4. Correlation coefficients for estimators of the peak-peak correlation function of peaks above la in the MDM simulation shown 
in Fig. §. 
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Figure 5. Mixed dark matter model correlation functions with errors expected from MAP (left) and Planck (right). The MAP simulations 
have a 12.6 arcminute beam, and Planck 5.5 arcminutes. Pixel sizes are taken to be 5 and 2 arcminutes. The simulations were done on 
square grids of side 3072 pixels, and the error bars (from run-to-run variations) rescaled to mimic 75% sky coverage. Note that the large 
pixels in MAP introduce discreteness errors for separations up to about 20 arcminutes, and the theoretical MAP curve is reliable only 
for separations greater than 15 arcminutes. 

required for the gaussian test is just the power spectrum; it does not rely on a satisfactory parameter estimation. For practical 
purposes, there are some issues still to address; detector noise can by easily incorporated into the analysis provided it is 
gaussian. Practical de-striping algorithms (e.g. Delabrouille 1998b, Revenu et al. 1998) to reduce 1/f noise, or methods to 
remove sidelobe contamination (Delabrouille 1998a) will need to be tested numerically to ensure that residual non-gaussian 
contributions are small. One interesting possibility for dealing with residual point sources in the maps is to put some constraint 
on the curvature at the peak (through the parameter x) to exclude sharply-peaked sources. The same trick could be used 
to exclude very flat-topped peaks whose precise locations may be difficult to determine accurately after dust foreground 
subtraction. Of course, both these techniques can be applied to the data as well as incorporated into the theory. There is an 
alternative to this procedure for point sources uncorrelated with the CMB; the peak distribution is then a superposition of 
two point processes, and provided the flux and clustering properties of the point sources are known, their contribution to the 
power spectrum can be subtracted. The gaussian test made on the joint distribution, as the correlation function will be a 
suitably- weighted combination of the two components. An aside which is probably not relevant for this application is that the 
analysis applies equally to any monotonic local functions of 2D gaussian fields, as nothing changes except the interpretation 
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of the threshold, although in general the non-gaussianity would be revealed in the one-point temperature distribution. In 
addition, one might expect any transformation process to be local in 3D, in which case the mixing of information through the 
last scattering surface would destroy the local nature. 



APPENDIX 

Joint distribution of variables 

Consider two points separated by r. For convenience, place them on the x-axis. In the spirit of Regos & Szalay (1995), we 
define complex quantities 



2/to( x ) = — / d 2 \i 5k exp(ik • x) exp(im(9)fc r: 



(1) 

where 8 is the angle between k and the a;-axis. With this definition, y" m (x) = y^*(x). Then we require the following complex 
variables: 

o X 
002/0 = o 

012/1 = fix + ifiy 

022/0 = fixx + Syy 

022/2 = fixX + Syy + 2l5 X y 



(2) 



The correlations of these (at the same point) are 



d 2 kd k'(<5k<5k') exp(ik ■ x — ik' • x) exp(im8 — im'9')k n ^' 



(3) 



With the definition of the power spectrum (hi), the integral simplifies to 



<2/m(x)s/m*(x)) = 



2 s-K 
T (n + n')/2°mm' 



(4) 



where 8 K is the Kronecker delta. The cross-correlations between variables at x and at x + r are: 



<2/m(x)t/w*(x + r)) 



d 2 kd 2 k'(5k<5k') exp(ik • x — ik' • (x + r)) exp(im8 — im'd')k n+ ' 



(5) 



Inserting the power spectrum removes the k' integration, leaving a 2D integration over k. The angular part of this is 



/■2-rr 

Ib = / d8 exp(— ikr cos 8) [cos(m — m')8 + i sin(m — m')#] 
Jo 



- 27ri |m_m,| 7i 

i " t ° \m — m' | 



(-Ax) 



(6) 



where J is a Bessel function, and J/(-z) = (-l)V;(z). Thus (y™ (x + r)y™ ,* (x)) = (-l) m + m ( y « ( x )y™ ,*(x + r)) and we find 
the cross-correlation matrix 

/ Aooo — A020 —Aon —Aon 

— A020 A220 A121 A121 

Aon -Am —A110 — A112 

Aon — A121 — A112 A110 

A022 — A222 — A121 — A123 

A022 — A222 — A222 A121 



(7) 



A022 A022 

— A222 — A222 
— A121 A123 

A123 — A121 
A220 A224 
A224 A220 / 

where the order of the variables is {j/q, y 2 , y\, y^ 1 , y\, y^ 2 } with quantities evaluated at x down the left, and at x + r along 
the top. 

(27T) 3 



dkk 1+a+0 P(k)J v (kr) 



We now make real variables 



(8) 



v = 2/o 
X = 



2 
</o 



Y = i(2/2 2 +2/- 2 ) 



© 0000 RAS, MNRAS 000, 000-000 



8 Alan F. Heavens and Ravi K. Sheth 



z 



1 I 2 2 s 

-^\V2 -y-2) 
2> 



0) 

These definitions agree with (0) in the main text. The last two are (85/ dx) jo\ and (OS / 'dy) / 'o\. The entire correlation matrix 
for the variables in the order (i>\,X\,r) x \, Y\, V2, X-z, Tj X 2, Ya,r] y i,r]y2, Z\, Z2) splits into an 8 x 8 matrix (top left) 



Ah 





f 1 


7 






Aooo 


Ao20 


— A011 


A022 


\ 




7 


1 






A020 


A220 


— A121 


A222 










1/2 




A011 


A121 


(Ano - 


All2)/2 


— (Al21 — Al23)/2 












1/2 


A022 


A222 


(A121 - 


Al 23 )/2 


(A220 + A224)/2 






Aooo 


A020 


Aon 


A022 


1 


7 












A020 


A220 


A121 


A222 
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1 












—Aon 


— A121 


(Ano — An2)/2 


(A121 — Ai23)/2 






1/2 








\ A022 


A222 


-(Al21 — Al23)/2 


(A220 + A224V2 










1/2 


/ 


a 4 x 4 block (bottom- 


right): 


















f 


1/2 


(Ano + An 2 )/2 




-(A 


121 + A 


123)/2 \, 










(Ano 


+ An 2 )/2 


1/2 


(A121 + Ai23)/2 


















(A121 + Ai23)/2 


1/2 


(A22O — A224)/2 










^ -(A121 


+ Ai 23 )/2 




(A22O — A224)/2 




1/2 


J 









(10) 



(11) 



Off-diagonal blocks are all zero. The 4x4 matrix is readily inverted. We used Mathematica to form the quadratic ^^(Mg)" 1 ^ 
(where v is the relevant part of the variable list, with the first derivatives set to zero). Since it produces several thousand 
lines of Fortran, it is not reproduced here. 
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